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Abstract 

In this work, we study the out-of-equilibrium many-body tunneling dynamics of a Bose-Einstein 
condensate in a two-dimensional radial double well. We investigate the impact of interparticle 
repulsion and compare the influence of angular momentum on the many-body tunneling dynamics. 
Accurate many-body dynamics are obtained by solving the full many-body Schrodinger equation. 
We demonstrate that macroscopic vortex states of definite total angular momentum indeed tunnel 
and that, even in the regime of weak repulsions, a many-body treatment is necessary to capture 
the correct tunneling dynamics. As a general rule, many-body effects set in at weaker interactions 
when the tunneling system carries angular momentum. 
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I. INTRODUCTION 


The experimental realization of Bose-Einstein condensates (BECs) in trapped dilute ul¬ 
tracold gases CHS] paved the way for studying the dynamics of many-boson systems. In par¬ 
ticular, tunneling phenomena with BECs have been of interest in recent years. A prominent 
example is the tunneling behavior in double-well potentials, often termed bosonic Josephson 
junctions in this respect [6j, which has been theoretically predicted um and observed in 
experiments [9* HQ]. Nowadays, the full many-body Schrodinger dynamics of a repulsive 
BEC in a one-dimensional double-well potential is available and shows the development of 
fragmentation, indicating that a many-body treatment is necessary in order to obtain the 
correct dynamics mm- 

Tunneling phenomena in two-dimensional (2D) trapped BECs, which is the main focus 
in this work, have been addressed recently. Especially the tunneling dynamics of trapped 
vortices were studied, e.g., in an harmonic potential with a Gaussian potential barrier [13], 
in 2D superfluids ra, between two Gaussian wells [15j, or between two pinning potentials 
eg. In three-dimensional double-well potentials, macroscopic superpositions of vortex states 
during the tunneling dynamics have been found m- 

The purpose of this work is however to investigate the full many-body Schrodinger dy¬ 
namics of a tunneling 2D BEC with definite total angular momentum. To this end, we 
consider a 2D radial double-well trap (see Fig. [Tj) . We discuss repulsive condensates made 
of N = 100 bosons with zero total angular momentum, L — 0, and vortex states with total 
angular momentum L = N. We demonstrate numerically that BECs carrying definite total 
angular momentum do indeed tunnel through the potential barrier. We compare the impact 
of angular momentum on the tunneling process and show that many-body effects set in 
at weaker interactions when the tunneling system carries angular momentum. A general 
conclusion stemming from our results is that the long time tunneling dynamics of 2D BECs 
cannot be described by a standard mean-field, like the Gross-Pitaevskii (GP) equation, even 
in the regime of weak interaction between the bosons. Thus, one is in need of an accurate 
many-body theory. Our tool of choice is the multiconfigurational time-dependent Hartree for 
bosons (MCTDHB) method [T8B2I] which has been applied [221 [27] and benchmarked [28] 
in various numerical studies on repulsive many-boson systems in recent years, in particular 
for BECs in 2D traps, see Refs. [29b35[ . 
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The paper is structured as follows. In Sec. |TT[ we present the theoretical framework of 
our study, i.e., we describe the system’s setup and introduce analysis quantities relevant for 


our purpose. The results are shown in Sec. Ill, where we split the discussion of the static 


considerations of the system (Sec. Ill A) from the tunneling dynamics of L = 0 (Sec. Ill B) 


and of vortex states (Sec. Ill C). Concluding remarks are given in Sec. IV We provide 
additional information on the numerical preparation of vortex states at the many-body level 
as well as a brief discussion on the numerical convergence in appendices [A] and [B] 


II. THEORETICAL FRAMEWORK 


A. Hamiltonian and setup 


The equation-of-motion for ultracold bosons in the 2D circular trap is the time-dependent 
Schrodinger equation 

i^ t m = H\*) : (i) 

where T) is the many-body wave-function which depends on the coordinates of all particles 
and time. The many-body Hamiltonian H is given by 

N N 

H = £ kn) + J2 tdR - >V|), (2) 

i =1 i<j=l 


with the single-particle Hamiltonian h(f) = — |A + V(r), r = (r, 6*), comprised of the 
kinetic energy and the external potential V, and the two-body interaction W. We consider 
dimensionless units that are obtained by dividing H by where h is Planck’s constant, 
d is a length scale, and m the boson mass. For typical realistic parameters see Ref. [fJTJ. 

The short-range interaction between the bosons is modeled by a Gaussian function [36, 


[37J, 


W{\fi — fj\) = 


- X ° .-|k-r,P/2 0 


27rcr 


(3) 


with cr = 0.25. To quantify the interaction strength, we introduce the mean-field interaction 
parameter A = Xo(N — 1) which is chosen to be positive in the following to describe repulsion 
between the particles. 

The external potential for the dynamics describes a 2D circular crater with a ring-shaped 
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central potential barrier, forming a 2D radial double well. It explicitly reads 


\ B e -2 ( r-i?s ) 4 + C e -0,5 b-«c) 4 if r < R c 

V(r) ={ (4) 

1C if r- > 


where Rb and Rc are the radial positions of the barrier and the crater’s wall, and B and 
C their heights. Throughout this work, we set B = 1 and the crater’s wall is kept fixed 
at Rc = 9.0 with constant height C = 200. A schematic plot of the setup is shown in 
Fig-El The chosen geometry serves as a natural way to combine tunneling phenomena in a 
double-well system with the conservation of angular momentum. Experimental realizations 
of this kind of trap have been achieved recently, see Refs. |541 53] . 

The time-dependent Schrodinger equation Eq. ([Tj) with the full many-body Hamiltonian 
Eq. ([2]) is solved by applying the MCTDHB(M) method. Its main idea is to express the 
wave-function of the system by a superposition of permanents {|n; t)} comprised of M time- 
adaptive orbitals {0j(r, t) : 1 < j < M}, 


W)) = | n;t), 

n 


(5) 


where n = (ni,..., um) is a vector carrying the individual occupation numbers of the orbitals 
and {Cfi(t)} are the expansion coefficients. It is important to note that both the expansion 
coefficients and the basis set are time-adaptive and determined by the Dirac-Frenkel, time- 
dependent variational principle. For M — 1, the MCTDHB theory boils down to the GP 
theory which is the standard mean-held commonly used to describe time-dependent BECs. 
Further details on the MCTDHB(M) method are given in the literature [T8ff2T| . We use the 
implementation in [38]. The simulations in this work are performed on a square box of size 
[-12,12) x [—12,12) with 128 x 128 grid points. A translation into real units can be found 
in Ref. [39]. The obtained results are converged to the accuracy given below. 
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B. Quantities of interest 


In order to analyze the properties of the system, some useful quantities are needed to be 
introduced. The reduced one-particle density matrix of a many-body system is given by 

M 

p (1) (r|r / ; t) = (T|T t (f / , f)T(r, f)|T) = p jk (t) t)(f) k (r, t) = 

j,k= 1 
M 

= n k (t)a* k (r ', t)a k (r, t), (6) 

k =1 

with the annihilation operator dh(fy t) = JT bj(t)(/)j(f,t ) that annihilates a boson at position 
r at time t. Its conjugate counterpart denotes the creation operator, creating a boson at po¬ 
sition r at time t. The eigenvalues {n k (t)} of Pj k (t) = Yhn' n t\oj(t)b k (t)\n] t) 

are called the natural occupations and the eigenfunctions the natural orbitals of 

the reduced one-particle density matrix. The natural occupations are arranged in descend¬ 
ing order, i.e., ni(t) > n 2 (t ) > ... > n M (i), with n j(t) = N. If there is only one 
macroscopic eigenvalue, the system is said to be condensed jj4TJ], whereas the system is said 
to be fragmented if two or more eigenvalues are macroscopic mm- 
The density of the evolving system is given by the diagonal of Eq. (§- 

p(r;t) = p m (f\r\t). (7) 

For the following studies on the tunneling dynamics, the trap is separated into an internal 
and an external part by setting the ring-shaped potential barrier at the position Rb■ We 
refer to these subsystems as the IN and OUT regions (see Fig. [Tj) , representing the trap’s 
center and the external rim, respectively. The occupation probabilities of the two parts are 
defined as 

p m = p(f]t)df (8) 

iV Jr<R B 

and 

PovTit) = 1 - P m (t) = i f p(f]t)df, (9) 

JR B <r<R c 

respectively (the density practically vanishes for r > Rc )■ We will among others use these 
probabilities to follow the dynamics in time. 
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FIG. 1. (Color online) Schematic plot of the system’s setup. Shown is a cut through the potential 
V ( r ) (shaded gray area) of the 2D radial double well, as well as typical density cuts of the ground- 
state GP orbitals 0™ and 4 >< a JT for both L = 0 and L = N. The number of bosons is N = 100. 
The IN and OUT regions of the trap are separated by the potential barrier at Rb (vertical dotted 
black line). For L = 0, the density maximum in the IN region is located exactly in the trap’s 
center, i.e., at r = 0 (dotted red curve). On the other hand, for L = N, one can clearly see 
the node of (dashed red). In the OUT region, for L = N (blue, double dotted) is 

pushed further to the crater’s wall than for L = 0 (solid blue). Horizontal lines denote typical 
(angular-momentum-dependent) values of the ground-state energies at the mean-field level for the 
corresponding orbitals, i.e. E™ (solid red) for </>™ and T ( s °lid blue) for ^>° UT . See text for 
more details. All quantities are dimensionless. 

III. RESULTS AND ANALYSIS 

A. Static considerations 

Before we start to study the dynamics, we want to analyze some static properties of the 
system which we will later use. Therefore, we first investigate the dependence of the ground- 
state energy on the position of the barrier Rb for N = 100 bosons at the GP, mean-held 


6 
















level (M = 1). We use the same methodology as in Ref. [3Tj for L — 0 and Ref. [35j for 
vortex states. The bosons are either trapped in the IN or OUT region separately. We would 
like to study the influence of the interaction strength A and the total angular momentum L. 
We call the corresponding energies E™ and E^ VT and denote the respective ground-state 
orbitals as and ( see Fig. [l] for more details). The respective ground states have 

been obtained by propagating in imaginary time, see Refs. HSUS] for further details. 


1. Impact of the interaction A 

The ground-state energies for different barrier locations Rb are depicted in Fig. [2^i. The 
total angular momentum is zero, L = 0. We consider the non-interacting single-particle 
case (A = 0) and the weakly-interacting GP case (A = 2). A discussion on the repulsion 
strengths in real units can be found in Ref. [56]. In both cases, there are certain radii where 
the ground-state energies for the IN and OUT regions coincide. These radii are termed the 
crossing points or simply the crossings in this work. They clearly depend on the repulsion 
strength. In the non-interacting case, we observe the crossing point at R\ = 3.271, whereas 
for A = 2 the crossing point is at a larger radius, R 2 = 3.412. These two radii are used in 
the dynamics below. 

We did the same simulation for ten times weaker interaction, A = 0.2, and three times 
stronger interaction, A = 6, and obtained, respectively, crossings at R B = 3.287 and R B = 
3.608 (not shown). We conclude from this that an increase of the repulsion strength A shifts 
the crossing point to a larger radius. From the viewpoint of trapping the particles in the IN 
region, this can be explained as follows. Due to the repulsion, the particles tend to separate 
from each other and the corresponding orbital is broadened. To account for the additional 
space needed in the trap’s center, the barrier needs to be shifted further towards Rc , he., 
to a larger radius. 


2. Impact of the angular momentum L 

In the previous subsection, we have shown that increasing the repulsion strength leads 
to larger values for the crossing point. The results were obtained for zero total angular 
momentum (L = 0). Here, we would like to study the impact of angular momentum L. 
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Details on the numerical preparation of such vortex states as well as on the conservation 
and measurement of L (during the time evolution) are given in Appendix [Aj 

Fig. ^ shows the ground-state energies where the total angular momentum of the system 
is L = N. The crossing point is located at R 3 = 4.089 (A = 0). This already shows that 
angular momentum strongly affects the position of the crossing point, comparatively stronger 
than the interaction strength. For a rough estimate, the crossing point for a relatively strong 
interaction A = 20 with L = 0 is at Rb ~ 3.95, which is still smaller than P 3 . Another 
observation is the significant increase in energy compared to the system with L — 0. The 
reason why angular momentum shifts the crossing point further towards the crater’s wall 
can be explained by considering the centrifugal barrier originating from the kinetic energy 
operator in 2D. It acts like an additional repulsive potential, pushing the particles away 
from the trap’s center. It also leads to the characteristic nodal structure of the respective 
GP orbital 0™ (see dashed red curve in Fig. [I]). 


B. Tunneling dynamics for L = 0 


1. The non-interacting system 


In this subsection, we analyze the tunneling dynamics for the non-interacting system 
(A = 0), i.e., the single-particle case where the total angular momentum is zero, L = 0. 


Thus, the orbital angular momentum is also zero, i.e., I = 0 [see Eq. (Al) in Appendix 
[A] for more details]. Fig. [3^ depicts the time evolution of Pout(^), Eq. ([9]), for different 
barrier positions Rb- The particles were prepared in the OUT region initially. In all cases 
considered, the bosons tunnel in a periodic manner between the two parts of the trap. The 
period and amplitude of the density oscillations clearly depend on Rb- Close to the crossing 
point Ri = 3.271 (see Fig. §>), the amplitude is maximal; almost the whole density is 
involved in the tunneling process. Also the period of a single tunneling cycle is maximal. If 
one leaves R ,\, both the amplitude and the period become smaller. 

Interestingly, the dynamics at the crossing point seem to be not sensitive to the initial 
condition, i.e., it does not matter whether the bosons are released from either the IN or the 
OUT region. This can be deduced from Fig. [3 ]d where the impact of the barrier location Rb 
on the period r and amplitude A of the P IN (t)- and Pout (^-oscillations is shown. Whereas 





FIG. 2. (Color online) (a): Dependence of the ground-state energy on the position of the radial 
barrier Rb ■ The angular momentum is L = 0, the number of bosons is N = 100 and the number 
of orbitals used is M = 1. E™ (dashed red curve) and Ay T (solid green) denote, respectively, 
the ground-state energies for the IN and OUT regions for A = 0, E™ (dotted blue) and Ag T 
(magenta, double dotted) for A = 2. The black vertical lines at R\ = 3.271 and R 2 = 3.412 denote 
the crossing points for A = 0 and A = 2, respectively. Increasing the repulsion strength shifts the 
crossing point to larger radii, (b): Same as in (a) but for L = N. The corresponding crossing 
points are R 3 = 4.089 (A = 0) and R 4 = 4.093 (A = 0.2). The angular momentum L affects the 
position of the crossing points comparatively stronger than the repulsion strength A. Note that the 
angular momentum is macroscopic; all particles carry it. See text for more details. All quantities 
are dimensionless. 
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for Ri there are no striking differences observable, the initial condition starts to matter as 
soon as the barrier location is varied. 

The period of the tunneling oscillations serves as a meaningful characteristic time scale 
for the dynamics. At the crossing point Ri, we obtain T\ = n/R = 110.231 with R = 
(0(? UT \M<t>™)\ = 2-85 • 10~ 2 . This is in a very good agreement with the result from Fig. [ 3 } 
R has a similar physical meaning as the hopping parameter in the Bose-Hubbard model, 
see, e.g., Ref. [12j. ft can be seen as a measure for the (unnormalized) transition probability 
amplitude from the state (0™) to the state |0° UT ) under the influence of h. Thus, the larger 
R is, the faster a particle tunnels through the barrier. We will use the same time scale for 
the case of interacting bosons in the next subsection. 


2. The interacting system 

We now turn to the impact of the interparticle repulsion A on the tunneling dynamics 
for L = 0, both at the GP, mean-held level and especially at the many-body level. We split 
the discussion of the results to two interaction strengths, A = 2 and A = 6, which we term 
in what follows weak and strong, respectively. We point out that the distinction between 
‘weak’ and ‘strong’ implies a classification with respect to the underlying physics. Whenever 
two modes faithfully describe the dynamics, we refer to the interaction as weak, whereas the 
interaction is termed strong whenever more modes are needed. Further physical distinctions 
between the dynamics with weak and strong interactions are discussed below. 

It will turn out that for the many-body dynamics, numerical convergence for the regime of 
weak repulsion is reached with M = 4 time-adaptive orbitals. We benchmark this result in 
Appendix [B] Due to the observations from the non-interacting system, we will focus on the 
tunneling dynamics at the corresponding crossing points. To recall, these radii depend on 
the interaction strength. 

To study the impact of weak repulsion, we set A = 2 and put the barrier at i? 2 = 3.412 
(see Fig. [2^). We obtain similar results for starting either from the IN or the OUT region 
and will therefore only concentrate on the latter. Fig. [I] shows the time evolution of Pout(^) 
at the mean-field (GP) and many-body [MCTDHB(4)j levels. For the first cycle, the two 
approaches give very similar results. However, differences set in afterwards. On the one 
hand, the mean-field dynamics show essentially unperturbed oscillations which are very 
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FIG. 3. (Color online) (a) Single-particle tunneling dynamics for different barrier positions Rb- 
The orbital angular momentum is l = 0. The dynamics are started from the OUT region. Close 
to the crossing point R\ = 3.271, the amplitude of Pout(£) is maximal and almost the whole 
density tunnels (red solid curve). The period is in a very good agreement with the predicted 
value T\ = 110.231. Leaving the crossing point, both the amplitude and period become smaller 
(Rb = 3.1, dotted blue; Rb = 3.5, dashed black), (b) Dependence of the amplitude A and period 
t (plotted on the same axis) on Rb for the same parameters as in (a). The indices IN and OUT 
denote the two different initial conditions considered. Close to the crossing point R±, both the 
amplitude and period become maximal and give the same result for the two initial conditions. See 
text for more details. All quantities are dimensionless. 
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similar to the non-interacting case. Only the period is slightly shorter due to the repulsion. 
On the other hand, the oscillations in the many-body case show damping of the amplitude. 
Pout(^) saturates after approximately 14 cycles with 52 particles in the external rim. This 
is very close to the equilibrium distribution of confining exactly 50 bosons in each part of the 
trap. Deviations from this can be explained with the accuracy of our numerical procedure to 
determine the crossing points. Thus, the observed behavior is similar to the density collapse 
in a bosonic Josephson junction, see, e.g., Refs. There, the density oscillations 

between the two wells are suppressed after some time. 

The evolution of the natural occupations show that the many-body dynamics are domi¬ 
nated by two natural orbitals, i.e., the system evolves from being essentially fully condensed 
[n 2 (0) = 0(1(T 3 ), n 3 (0) and 774 ( 0 ) are even smaller] to being two-fold fragmented. After the 
density oscillations have collapsed, the occupation numbers of the first two natural orbitals, 
aq and ct 2 , are approximately 56.7% and 41.1%, respectively. The remaining two orbitals are 
significantly less occupied (approximately 1%). It is interesting to note that the curves for 
ni(t) and n 2 (t) tightly envelope the oscillations of Pout(^)> indicating that fragmentation 
and the collapse of the density oscillations are closely related, see Fig. [4] 

Let us have a closer look into the structure of the first two natural orbitals aq and ct 2 at 
t = 18.4Ti, i.e., when the density has already collapsed at the many-body level. The naive 
notion would be to have one localized orbital in the IN region and another one localized 
in the OUT region. However, Fig. [5] shows that they are rather delocalized in space, i.e., 
distributed over both subsystems. In contrast to that, the GP orbital is localized in the 
external rim and does not account for the occupation of the trap’s center. This indicates 
a major difference between the fragmented many-body system and the condensed system 
when described at the GP level. 

The above simulations have been repeated for ten times weaker interaction, A = 0.2, 
leading to qualitatively similar results. The only crucial difference at the many-body level is 
that it takes much longer (more than 200 cycles) until the density oscillations are completely 
suppressed. 

We now turn to the case of strong repulsion (A = 6) with crossing point Rb = 3.608. 
The occupation numbers at t — 0 are: 77i(0) = 99.37%, n 2 (0) = 77 . 3 ( 0 ) = 0.28% and 774 ( 0 ) = 
0.06%, i.e., the system is still essentially fully condensed. In the many-body simulation, we 
observe that during the time evolution all 4 natural orbitals become significantly occupied, 


12 


already after 4 cycles. We can therefore not claim that we reached numerical convergence 
with M = 4 time-adaptive orbitals in this interaction regime, and we will therefore not 
physically interpret the underlying natural orbitals as we did in Fig. [5]for weak repulsion. To 
reach numerical convergence, one would need to allow for additional time-adaptive orbitals, 
ffowever, already one additional orbital for the same number of particles would substantially 
increase the configuration space of size ( N+ such that the computational effort would 
exceed the scope of this work. Nevertheless, we can draw some conclusions. Starting as 
above from the barrier located at the crossing point, the time evolution of Pout(^) at the 
mean-held and many-body levels are depicted in the inset of Fig. |4} At first, we see that 
much less particles are involved in the tunneling process compared to the case of A = 2. 
Secondly, the periods of the oscillations for the mean-held and many-body descriptions are 
different. The equivalence of both approaches is broken already during the hrst cycle. 

Another interesting observation is that the choice of the initial condition becomes impor¬ 
tant for stronger repulsion, meaning that the dynamics of Pm{t) starting in the IN region 
and of Pout(£) starting in the OUT region are no longer equivalent (not shown). It is worth 
mentioning that this is very different compared to the (one-dimensional) symmetric double¬ 
well system. In the latter, one would never observe any differences between (t ), when 
starting in the left well, and Pii(t), when starting in the right well, no matter how large the 
interaction strength is. This is due to the perfect spatial symmetry of the external potential 
which is not given in the 2D circular trap. Since the two IN and OUT subsystems have 
different topologies and sizes, one cannot expect their level structures to be equivalent; This 
leads to different dynamics for the initial conditions because stronger interaction means that 
more than just the lowest-in-energy time-adaptive orbitals are involved. 


C. Tunneling dynamics of vortex states 

In the previous subsection, we have investigated the tunneling and fragmentation dynam¬ 
ics for bosons with zero total angular momentum, L = 0. Here, we report on the tunneling 
dynamics of vortex states with definite total angular momentum L = N, and explore both 
similarities and differences compared to the dynamics of L — 0. We will again discuss 
non-interacting and interacting bosons separately. 
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FIG. 4. (Color online) Mean-field (M = 1) and many-body (M = 4) tunneling dynamics of 
the interacting system at the crossing point R -2 = 3.412. The angular momentum is L = 0, the 
number of particles is IV = 100, and the interaction strength is A = 2. The bosons are initially 
prepared in the OUT region. Whereas the GP theory predicts unperturbed oscillations of -PoutW 
(dotted gray curve), the amplitude is damped at the many-body level (solid red) and saturates 
after approximately 14 cycles. Only the first two natural orbitals are significantly occupied (solid 
green and blue), the other two (solid magenta and light blue; curves lie atop of each other) carry 
just a small fraction of the particles. However, the frequencies of the GP and many-body density 
oscillations are essentially the same. Inset: Same simulation but for strong repulsion (A = 6) 
at the respective crossing point Rb = 3.608. Many fewer particles tunnel back and forth; the 
dynamics obtained from the GP and MCTDHB(4) descriptions start to deviate from each other 
already during the first cycle, both in their frequencies and amplitudes. After 4 cycles, the four 
natural orbitals are already macroscopically occupied (not shown). Note the different time scales. 
See text for more details. All quantities are dimensionless. 


1. The non-interacting system 


In this subsection, as a starting point, we investigate the dynamics of vortex states 
where the constituent particles do not interact, i.e., the single-particle case with orbital 
angular momentum l = 1 [see Eq. (Al) in Appendix [A] for more details]. We first study the 
impact of the barrier location Rb on the tunneling dynamics. We compare the periods and 
amplitudes of the density oscillations between the IN and OUT regions for several values of 
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FIG. 5. (Color online) Real, imaginary, and absolute value of the GP, mean-field orbital (top 
panels) and of the first two natural orbitals a\ and a - 2 of the many-body simulation (M = 4, 
middle and bottom panels) at t = 18.4 t±. The other parameters of the system are the same as in 
the main figure of Fig. |4j Whereas the mean-held orbital is localized in the external rim, a\ and 

from the many-body computation are delocalized, covering both the IN and OUT regions. See 
text for more details. All quantities are dimensionless. 

R b in Fig. i We recall that for l — 1 the crossing point in the non-interacting system is 
located at i? 3 = 4.089, which is larger compared to the non-interacting system with l = 0 
where the crossing point is located at R\ = 3.271 (see Fig. [2|. Now, for l = 1 and barrier 
position i?3, both the amplitude A and period r are maximal with values of A = 48 particles 
and r = 66.4. 

In order to calculate the characteristic tunneling time 72 of a vortex state, we prepared 
the two ground-state orbitals with l = 1 for the IN region (0 q N ) and for the OUT region 
(0o UT , cp. Fig. [l|, respectively, and calculated J 2 = |(0o UT |^|0o N )l = 4-704 • 10~ 2 . This in 
turn gives t 2 = 7 r/J 2 = 66.786 which is in a very good agreement with the result from Fig. [6j 
Compared to the case of l = 0, the vortex state tunnels almost twice as fast between the 
two parts of the trap. This can be anticipated from its higher energy and the lower effective 
barrier under which it has to tunnel through. We recall that the barrier height is kept fixed 
at B = 1 . Moreover, the dynamics give similar results at R 3 for the two initial conditions 
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FIG. 6. (Color online) Dependence of the amplitude A and period r (plotted on the same axis) on 
Rb in the non-interacting system, i.e., the single-particle case. The orbital angular momentum is 
1 = 1. The indices IN and OUT denote the two different initial conditions considered. Close to the 
crossing point R 3 = 4.089, the period and amplitude of the density oscillations are maximal and 
the two initial conditions give equivalent results for the tunneling dynamics. Leaving the crossing 
point, both the amplitude and period become smaller. See text for more details. All quantities are 
dimensionless. 

considered. Our study of the interacting vortex states will therefore focus on the dynamics 
at the corresponding crossing points where the bosons are initially prepared in the OUT 
region. 


2. The interacting system 

We now investigate the impact of the interparticle repulsion A on the tunneling dynamics 
of a vortex state with L = N. It turns out that at the many-body level, vortex states are 
more sensitive to the strength of the repulsion than states with L = 0 (see details below). In 
other words, many-body effects set in at weaker interactions. Thus, we split, the discussion 
for the regime of weak repulsion (A = 0.2) and strong repulsion (A = 2). Note that the 
regimes for L = N are different than for L = 0, see the discussion in Sec. |IV| 

For A = 0.2, we put the barrier at R 4 = 4.093 (see Fig. [2 |d) and prepare the bosons in 
the OUT region. At the mean-held level, the vortex state tunnels essentially unperturbed 
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between the IN and OUT regions, i.e., there is no damping of the amplitude observable. 
The result is thus very similar to the one from the non-interacting case and is therefore not 
discussed further (see upper inset of Fig. [7]). 

The many-body dynamics on the contrary are more intriguing. Fig. [7] shows the MCT- 
DHB(4) results for the time evolution of Pout(^) together with the corresponding natural 
occupation numbers. The initially-coherent system [ 712 ( 0 ) = (9(10~ 8 )] evolves to a two-fold 
fragmented condensate where the two time-dependent natural orbitals both carry angular 
momentum 1 = 1. The remaining two orbitals do essentially not participate in the dynamics 
since their occupations stay below 0.1% throughout the time evolution. 

The evolution of Pout(^) shows damping of the oscillations. The damping resembles the 
case of L = 0 where we observed the same phenomenon. The behavior at longer times is 
such that the system tends to distribute the particles equally between the IN and OUT 
regions, and the density oscillations get suppressed. The GP description does not predict 
the density collapse of a vortex state and thus a many-body treatment of the dynamics is 
necessary. This conclusion is supported by Fig. [8] which shows the structure of the first two 
natural orbitals from the many-body computation as well as the corresponding mean-field 
orbital at an intermediate time step t ~ 183.3 r 2 - The GP theory does not account for the 
occupation of the trap center and only predicts particles in the external rim. I 11 contrast to 
that, the delocalized orbitals of the many-body computation show significant occupations 
of both the IN and OUT regions. In the trap’s center, the characteristic nodal structure of 
vortex states can be observed. The second orbital a 2 has an additional ring-shaped node, 
matching the location of the barrier at f? 4 . 

We repeated the same simulation for strong repulsion A = 2 at the respective crossing 
point Rb = 4.123. By observing that already for this interaction strength all 4 natural or¬ 
bitals become macroscopically occupied, we deduce that for this value of A one has already 
left the weak repulsion regime for vortex states. Thus, vortex states are much more sensitive 
to interparticle repulsion in comparison to states with L = 0. We again stress that it is nec¬ 
essary to include more time-adaptive orbitals to reach numerical convergence. Nevertheless, 
already M = 4 orbitals are enough to expose substantial differences between the mean-held 
and many-body predictions already after a few cycles, both concerning the amplitude and 
frequency of the density oscillations (see lower inset of Fig. [7J note the different time scales). 
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FIG. 7. (Color online) Many-body (M = 4) tunneling dynamics of a vortex state at the crossing 
point = 4.093. The number of particles is N = 100, the total angular momentum is L = N, and 
the interaction strength is A = 0.2. The bosons are released from the OUT region. The damped 
oscillation of -Pout(£) (solid red) is enveloped by the natural occupations of a\ and a 2 (solid 
green, blue). The remaining natural orbitals do not become significantly occupied (solid magenta, 
light blue; curves lie atop of each other). Upper inset: Dynamics of Pout(^) for the GP, mean- 
field (M = 1, dotted gray) and many-body (M = 4, solid red) descriptions between t = 150 t\ and 
t = 155 Ti. The system’s parameters are the same as in the main figure. The mean-field description 
does not account for the density collapse; The vortex state is already fragmented. However, the 
tunneling frequencies of the GP and MCTDHB(4) results are essentially the same. Lower inset: 
Same simulation for strong repulsion, A = 2, at the respective crossing point Rb = 4.123. The 
other parameters of the system are the same as in the main figure. The dynamics obtained from 
the mean-held (dotted gray) and many-body descriptions (solid red) start to deviate from each 
other already after a few cycles, both in frequency and amplitude. After 10 cycles, the four natural 
orbitals are already macroscopically occupied (not shown). See text for more details. All quantities 
are dimensionless. 

IV. CONCLUDING REMARKS 

In conclusion, in the present work we have studied static properties and particularly the 
out-of-equilibrium tunneling dynamics of BECs and vortex states in a 2D radial double well. 

On the statics side, we showed that angular momentum and repulsion between the bosons 
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FIG. 8 . (Color online) Real, imaginary, and absolute value of the GP, mean-field orbital (top 
panels) and of the first two natural orbitals a\ and a -2 of the many-body simulation (M = 4, 
middle and bottom panels) at t = 183.3 T 2 - The system’s parameters are the same as in the main 
figure of Fig. [7j Whereas the mean-held orbital is localized in the external rim, the natural orbitals 
in the many-body computation are delocalized, covering both the IN and OUT regions. ot\ only 
has a node in the trap’s center, whereas «2 has a second, ring-shaped node at the barrier’s position 
R 4 . See for comparison Fig. [5] and the text for more details. All quantities are dimensionless. 

affect the location of the ground state, i.e., whether it is energetically favorable for the bosons 
to occupy either the IN or the OUT region of the trap. For a certain critical radius, which 
we termed the crossing point, the ground-state energies of the IN and OUT subsystems are 
equivalent. The position of the crossing point depends on both angular momentum and 
repulsion strength. The impact of angular momentum is however significantly stronger, i.e., 
it shifts the crossing point to much larger radii. 

On the dynamics side, we observed several similarities as well as clear differences between 
the tunneling of BECs with L — 0 and of vortex states with L = N. At first, we demon¬ 
strated periodic, Rabi-like tunneling between the IN and OUT regions in the non-interacting 
and weakly-interacting regimes for both states with L — 0 and vortex states with L = N. 
The density oscillations are most pronounced when the barrier is located at the correspond¬ 
ing crossing points. Leaving the crossing points leads to density oscillations with shorter 
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periods and smaller amplitudes, i.e., less particles are involved in the tunneling process. 
However, by comparing the characteristic tunneling times, we found that the tunneling of 
a vortex state with L = N, which is higher in energy, takes place on a (much) shorter time 
scale, being almost twice as fast as for states with L — 0. 

For both values of angular momentum considered, the time-dependent GP equation fails 
to describe the long-time tunneling dynamics as soon as the particles interact. For weak 
repulsion, both systems with L = 0 and L = N become essentially two-fold fragmented. 

The development of fragmentation is accompanied by damping of the density oscillations 
between the IN and OUT regions. The stronger the interaction A is, the faster the density 
oscillations are suppressed. The regime of weak repulsion for L = 0 is apparently more 
extended, which in turn means that vortex states are more sensitive to repulsion. We deduce 
this from the fact that already for A = 2 (which corresponds to weak repulsion for states with 
L = 0) more than two natural orbitals are significantly occupied in the many-body dynamics 
of vortex states with L = N at long simulation times. We have additionally shown that 
for the strongly repulsive systems (A = 6 for L = 0 and A = 2 for L = N ), at least 4 
time-adaptive interaction-dressed orbitals in the MCTDHB theory are necessary in order 
to describe the many-body tunneling dynamics faithfully. In both cases, the many-body 
dynamics quickly start to deviate from the corresponding GP predictions. 

The present work shows that the tunneling dynamics of BECs and of vortex states in a 
2D radial double well is many-body in nature. Furthermore, many-body effects set in at 
even weaker interactions when the tunneling system carries angular momentum. Obviously, 
with increasing interaction more and more many-body excited states are involved, where 
both radial and angular excitations can combine to assemble states of definite total angu¬ 
lar momentum. These may allow for an even richer and more intricate out-of-equilibrium 
dynamics than reported here. 
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Appendix A: Further details on vortex states and total angular momentum 

In this appendix, we give additional details on the numerical preparation of vortex states 
with L = N, as well as on the conservation of the angular momentum during the propagation 
in time of the many-body wave-functions. 

The numerical preparation of a vortex state in the 2D radial double-well potential can be 
achieved by multiplying radially symmetric, real-valued functions fj(r,t), 1 < j < M, with 
a phase of integer value: 

t) = fj (r) e llj6 , 1 <j< M, (Al) 

where 6 is the phase and lj the integer angular momentum per particle. The orbitals 
{(pj{f,t)} are then propagated in imaginary time, leading to the initial vortex states from 
which we study the dynamics in the main text. In principle, the individual orbitals involved 
can have different values lj. In this work, we start from essentially condensed systems and 
we thus concentrate on states where the underlying orbitals have the same lj = / = 1. 

During the time evolution, we measure the angular momentum of the underlying natural 
orbitals via 

4 = (&k(t)\L z \a k (t)), 1 < k < M. (A2) 

The angular momentum operator in z-direction is defined as L z = —i (x-§^ — For the 

regime of weak repulsion, we observe that Ik = 1 for the macroscopically-occupied orbitals 
is nicely conserved throughout the whole time evolution. The total angular momentum L is 
measured in our computations via the quantity 

L ~ y ' (Lz)jk(t) Pjk(t)i (A3) 

j,k 

with ( L z )jk(t ) = (<j)j(t)\L z \<j)k(t)). L serves as a good quantum number for the dynamics of 
both values of angular momentum considered in the main text. 

Our study suggests the implementation of a projection operator onto good-total-angular- 
momentum-eigenstates of the full many-body Hamiltonian Eq. |2| as a useful numerical 
development, extending our current method of preparing vortex states. In such a way, one 
would be able to prepare the system in many-body eigenstates of, e.g., non-integer total 
angular momentum per particle, L/N. This would possibly allow for a more sophisticated 
study of vortex states’ dynamics, in particular for stronger interactions, as well as the study 
of the evolution in time of the spatially-partitioned many-body vortices proposed in [35] , 


21 


Appendix B: Numerical convergence for L = 0 


The purpose of this appendix is to justify the restriction to M = 4 time-adaptive orbitals 
within the MCTDHB(M) theory in case of weak repulsion for L — 0. 

Fig. [9] shows the tunneling dynamics for N — 10 particles, allowing for M = 10 time- 
adaptive orbitals. The interaction parameter has been set to A = Ao (N — 1) = 2 as in the 
main text (see Fig. [4]). The tunneling dynamics are again dominated by two natural orbitals 
which carry the majority of the occupation probability. In addition to that, there are two 
more natural orbitals, occupied by roughly 10% of the bosons. However, the most important 
observation is that the remaining 6 natural orbitals do not become significantly occupied, 
i.e., •rij > ^{t) < 0.15%. This proves that numerical convergence of the many-body tunneling 
dynamics in the 2D radial double well with M = 4 time-adaptive orbitals is reached. 
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